Harnessing adaptive bistable stiffness of hair-cell-bundle structure for broadband vibration applications

This study presents an initial study on the adaptive bistable stiffness of the hair cell bundle structure in a frog cochlea, and aims to harness its bistable nonlinearity that features a negative stiffness region for broadband vibration applications such as vibration-based energy harvesters. To this end, the mathematical model for describing the bistable stiffness is first formulated based on the modeling concept of piecewise type nonlinearities. The harmonic balance method was then employed to examine the nonlinear responses of bistable oscillator, mimicking hair cells bundle structure under the frequency sweeping condition, and their dynamic behaviors induced by bistable stiffness characteristics are projected on phase diagrams, and Poincare maps concerning the bifurcation. In particular, the bifurcation mapping at the super- and sub-harmonic regimes provides a better perspective to examine the nonlinear motions which occur in the biomimetic system. The use of bistable stiffness characteristics of hair cell bundle structure in frog cochlea thus offers physical insights to harness the adaptive bistable stiffness for metamaterial-like potential engineering structures such as vibration-based energy harvester, and isolator etc.

The cochlea present inside the inner ear is one of the primary auditory organs in which sound is transduced from acoustic energy into an electrical signal. Its sensory receptors are called hair cells and feature densely bundled hair structures 1 . The primary function of the hair bundle structure is to send biologically induced electrical impulse signals to the brain in response to the vertical oscillation produced by the traveling wave propagation on the basilar membrane of the cochlea, as shown in Fig. 1a. Typically, the hair bundle structures of the cochlear outer hair cells in auditory systems consist of multiple tiny long cylinders (e.g., approximately 100 in an ear) called Stereocilia that lean on each other with tip links in the hair cell, as shown in Fig. 1b. The primary function of hair bundles originates from the dynamics of the tip link (elastic gating spring) connected to transduction channel gates, as shown in Fig. 1b. Apart from the restoring force produced by the bending stiffness, the tip links in the hair bundles deliver an external force to the transduction channel 2 . Furthermore, calcium ion (Ca2 + ) concentration plays an essential role in impulse signal transduction because the channel's open probability is determined by the hair bundle displacement and calcium ion concentration. As illustrated in Fig. 1b, individual hair bundles protruding from the bottom surface of the hair cells oscillate within the fluid-filled cochlea. Typically, in mammals, including humans, one row of the inner hair cells is aligned along the length of the cochlea in parallel to three rows of the outer hair cells. Stereocilia of the inner hair cells are linearly arranged, whereas outer hair bundles are arranged in a V-shaped pattern. These different morphologies are likely to reflect the distinct functions of inner and outer hair cells. Outer hair cells can amplify small oscillations, and therefore, significantly enhance the sensitivity and dynamic range of hearing. In contrast, inner hair cells do not amplify but transmit electrical signals to the auditory-nerve fibers 3,4 .
As shown in Fig. 1c, negative stiffness is experimentally observed when measuring the mechanical properties of the hair bundle structures of a frog [5][6][7] . The dynamic oscillator with this nonlinearity is referred to as a bistable oscillator, and it has a unique double-well restoring force potential and provides three distinct operating regions depending on the input amplitude. This bistable oscillator may be excited to exhibit aperiodic or nonlinear vibrations between the two wells (i.e., the negative-stiffness region). Thus, numerous studies have attempted to harness the bistability featuring a snap-through (also known as buckling) action in various engineering applications, such as energy harvesting devices, acoustic transducers, and fluidic sensors [8][9][10][11] .
To date, it has been hypothesized that bistability is harnessed to amplify mechanical auditory stimuli in hair bundle structures. The sensitivity of hearing can be enhanced by the mechanoelectrical phenomenon observed in the outer hair bundle structures, which can amplify hair cell motion over broad frequency ranges. This high sensitivity is also hypothesized to exist because of the combination of an adaptation and a negative stiffness property inherent in the hair bundle structures 12 . Hence, even though a sinusoidal stimulus is applied in a specific frequency range, the transferred sound signals can be amplified owing to the mechanical adaptation ability that shifts the region of highest sensitivity toward the active operation range of the hair cell structure. Although the detailed mechanism by which mechanical energy is amplified by the nonlinearity of hair cells describes well 3,5 , some relevant studies have still attempted to further identify the amplification mechanism [13][14][15] . In this study, we mimic only its mechanical property featuring a negative stiffness region and adaptation ability. The primary objective of this study is thus to investigate the nonlinear dynamic behavior of a bistable oscillator with adaptive bistable stiffness, mimicking hair cell bundle structure in a frog cochlea, and aims to harness its bistable nonlinearity that features a negative stiffness region for broadband vibration applications.

Mathematical formulation of hair cell bundle structure
In this study, the nonlinear stiffness characteristic of the hair cell bundle structure is represented by the nonlinear spring element of a 1-DOF vibrating system (i.e., a bistable oscillator), as shown in Fig. 2a. The nonlinear stiffness function of the hair cell bundle structure includes two identical positive and one negative stiffness coefficient. The equation of motion for a 1-DOF vibrating system can then be formulated as follows where F S (x) and F b (t) denote nonlinear spring and sinusoidal excitation forces, respectively. Similar to the stiffness curves displayed in Fig. 1c, the bistable nonlinearity is mathematically formulated using the hyperbolic tangent function to avoid instability due to the discontinuity of using the piecewise linear function. The effectiveness of the hyper tangent function for the sigmoid (S-shaped) function has been proved in many prior research 16 .
where the employed variables and symbols are as follows: k b(i) , such as k b1 and k b2 , stiffness values for each stage; x bpr = x b − x pr ; x pr , preload location; F S , total spring force; F S1 , spring force without the preload; F Spr , preload; F bp(i) , (2d) F bn(i) = x bpr + ν n(i) tanh σ b x bpr + ν n(i) − 1 . www.nature.com/scientificreports/ spring force on the positive side; F bn(i) , spring force on the negative side; ν p(i) , transition displacement on the positive side; ν n(i) , transition displacement on the negative side; σ b , smoothing factor. Thus, arbitrary piecewise-type nonlinearities of the bistable stiffness profiles can be determined by employing their numerical values. Case 2 in Fig. 2b shows the symmetric-type bistable stiffness, calculated by substituting the relevant properties. Because the stiffness of the hair cell bundle structure in Fig. 1c is observed by microscale [Force ( 10 −12 )/Displacement ( 10 −9 )] and it is an extremely small value ill-suited for real engineering applications, it was scaled to the macro level (i.e., 10 3 ), as listed in Table 1 11 . For different piecewise-type nonlinearities, multiple profiles can be produced using Eq. (2), as shown in Fig. 2b. In addition, to formulate each bistable stiffness curve, the smoothing factor σ b is set as 5 × 10 3 . Table 2 lists the relevant properties to determine three cases: Cases 1, 2, and 3. For example, Case 2 reflects the symmetric characteristics; however, Cases 1 and 3, which mimic the adaptive bistable stiffness, are shifted from the origin and become asymmetric. Thus, all possible bistable stiffness profiles can be defined by employing the proper values for F Spr and x pr based on the symmetric case, such as in Case 2. For example, Table 2 lists the values for F Spr and x pr to determine three different cases of nonlinear stiffness profiles, as shown in Fig. 2b. The proof mass m = 0.015 kg; damping coefficient c = 51.6 N s m −1 , assuming that the employed modal damping ratio, ζ, and natural frequency, ω n (f n ), are 5% and 516.4 rad s −1 (82.2 Hz) respectively. Additionally, the natural frequency ω n is obtained using the positive stiffness value k b2 = 4 × 10 3 N m −1 .

Nonlinear frequency response analysis bistable oscillator
To investigate the system nonlinear responses, the harmonic balanced method (HBM) was implemented in this study based on the Galerkin schemes because all the information of the system responses including the stability conditions in both time and frequency domains under the steady state conditions is efficiently obtained by employing the HBM and it has been known as one of the powerful tools to analyze the strongly nonlinear stiffness system [16][17][18][19] . First, the system response x b (t) and the input F b (t) can be considered as follows:  Transition displacement on the positive side ( www.nature.com/scientificreports/ where x m , x ak , and x bk are the mean and alternating parts of the cosine and sine functions for the system responses, respectively; F m is the average force (0.1 N); ω p and ϕ pk are the excitation frequency and phase angle (in this study, 0), respectively; η and k are the sub-and super-harmonic indices, respectively; N max is the maximum number of harmonics correlated with the harmonic index of the HBM. Fpk is the magnitude of the sinusoidal input force and selected to be 7 N such that it can induce the chaotic interwell vibrations 8 ; Assuming that the system is in a steady state, the Galerkin scheme in Eq. (1) is expressed as follows 20,21 .
Here, its relevant terms are defined as follows.
Furthermore, its nonlinear and input functions are defined as follows.
The relevant variables used are: ̟ t = ψ and ̟ = ω ω n , the non-dimensionalized time scale, and the normalized frequency value with the natural frequency ω n ; T = ητ , is the concerned period with respect to 0 ≤ t < T → 0 ≤ ψ < 2π ω n ; η is a sub-harmonic index; τ is the period of the fundamental excitation; k and l represent the incremental index where k = ω n , 2ω n , 3ω n · · · and l = 1, 2, 3 · · · . By employing the relationship between ẋ(t) = dx dt = ̟ dx dψ = ̟ x ′ and ẍ(t) = ̟ 2 x ′′ , the overall Galerkin scheme for the basic equation of Eq. (4) is expressed as follows: To determine the solutions of x bc and ̟ for each step, the Newton-Raphson method was implemented under the condition → 0 , where is considered a function of x bc and ̟ , such as x bc , ̟ . Detailed information on the derivation and descriptions of the HBM can be found in previous studies 16 .

Results and discussion
The HBM results with the root mean square (RMS) values of the displacement along with three different bistable stiffness curves are compared in Fig. 3, in which Cases 1 and 3 are asymmetric profiles, and Case 2 is a symmetric curve (a baseline for this study). Here, Cases 1 and 3 are shifted from the symmetric profile to the lower left and upper right directions, respectively, which corresponds to the mechanical adaptation capability. Cases 1 and 3 exhibit higher resonance values in the super-harmonic regimes than in the resonant regime of Case 2, as shown in Fig. 3b. However, Case 2 is highly affected by the sub-harmonic resonances compared with Cases 1 and 3, as shown in Fig. 3c. In addition, Cases 1 and 3 exhibit almost the same dynamic characteristics, whereas Case 2 F pk cos kω p t + ϕ pk . www.nature.com/scientificreports/ shows significant differences in the super-and sub-harmonic areas, except for the primary resonance. To obtain the nonlinear responses in Fig. 3, the employed values of η and N max for the HBM are 2 and 12, respectively. Figure 4 shows a comparison of the numerical simulation (NS) and the HBM results. In this study, the modified Runge-Kutta method was employed to obtain the NS result 22 . To reveal the sub-harmonic responses in greater detail, the values for η and N max of the HBM were 6 and 12, respectively. A greater number of complex sub-harmonic responses could be obtained with an increase in the number of sub-harmonic indices η , as shown in Fig. 4. The red dotted circles indicated as (A) and (B) represent the super-and sub-harmonic regimes, respectively, and the stability conditions were determined using Hill's method [23][24][25] . The details of the super-and sub-harmonic responses are shown in Fig. 4b and c. When the NS and HBM results were compared, the stable responses of the HBM correlated well with the NS results. However, the unstable response of the HBM was not correlated well because the unstable response is related to complicated system behaviors, such as quasi-periodic and chaotic phenomena. In addition, the differences between the NS and HBM results were due to different analytical processes. For example, NS is calculated based on the time domain integration using the initial conditions renewed at each prior and current step 22 . Thus, the NS can reflect all possible dynamic behaviors with time variations. However, the HBM determines its solutions efficiently by estimating the frequency and time domain information with the integer-based Fourier expansion, even though it cannot include all the possible time histories, such as transient responses [19][20][21]23,24 . Based on two different analytical approaches, the super-and sub-harmonic responses induced by bistable nonlinearities can be investigated in detail.
Nonlinear dynamic behaviors in the super-harmonic regime. The nonlinear dynamic characteristics of the NS and HBM were compared with the bifurcations in the super-harmonic regime, as shown in Fig. 5. In general, the bifurcation is defined as the dramatic changes of the system responses under the small variation of system parameters. For example, ̟ is considered as a parameter for the current system shown in Fig. 5. While the value ̟ is changed from the lower to the higher ranges, some of areas around ̟ = 0.3 show various solutions whenever each period is complete, which appears as scattered blue dots, as shown in Fig. 5. To calculate the bifurcations, the solutions of NS for each excitation frequency were obtained at the same orbital locations www.nature.com/scientificreports/  www.nature.com/scientificreports/ during 100 cycles of periodic motions after the transient responses were completely removed. When the stable solutions of the HBM are compared with those of the NS, they correlate well with each other, as shown in Fig. 5. However, unstable solutions of the HBM are closely related to bifurcation phenomena. By focusing on the area where the stability conditions change from unstable to stable (UTS), the UTS conditions reflect the bifurcation conditions well 18 . To analyze the dynamic responses in detail, the time histories can be reviewed in terms of two specific locations, ̟ = 0.3 and ̟ = 0.4, represented by (1) and (2), as shown in Fig. 5. For example, the system response at ̟ = 0.3 shows the unstable conditions from the HBM, corresponding to the period constituting the bifurcation cascade obtained by NS. In this regime, the time histories exhibit significant complexities, as shown in Fig. 6a. When the time histories calculated by both the NS and HBM are compared, the NS results include several harmonic terms rather than the HBM because HBM is constructed using integer-based increment, as described earlier. In addition, the harmonic components in terms of magnitude are clearly observed in Fig. 6c. While the FFT results from the NS, as shown in Fig. 6c, demonstrate multiple harmonic terms, the HBM shows only fundamental and 2nd harmonic terms. Meanwhile, the system responses at ̟ = 0.4 show good correlations between the two results between NS and HBM, as shown in Fig. 6b, because the dynamic behavior at this location is stable, which is also observed in the FFT spectrum shown in Fig. 6d.
In addition, the complexities of the dynamic behaviors could be efficiently examined from the phase diagrams and Poincare maps. Figure 7 shows a comparison of the phase diagrams and Poincare maps at ̟ = 0.3 and ̟ = 0.4. For instance, the phase diagram at ̟ = 0.3 includes more complex dynamic tracks, as shown in Fig. 7a. However, the phase diagram at ̟ = 0.4 shows only one clear cycle. The Poincare map in Fig. 7c demonstrates scattered points. However, the Poincare map at ̟ = 0.4 is concentrated at only one point.

Nonlinear dynamic behaviors in the sub-harmonic regime. The nonlinear dynamic behaviors in
the sub-harmonic regime are shown in Fig. 8. The dynamic behaviors in region (1) are affected by more nonlinearities because the stability conditions from HBM are determined to be unstable, and this region demonstrates severe bifurcation effects. However, region (2) shows relatively less complexity because this region pertains to the stable dynamic conditions determined by the HBM, even though the system responses in this region show a period-doubling effect. As shown in Fig. 9, the dynamic responses from the NS and HBM are compared based on the time histories. For example, Fig. 9a compares the time histories at ̟ = 1.1, marked as (1), in which the results from both NS and HBM are not well correlated. This is the same reason as that observed in superharmonic regions, implying that the system is affected by higher nonlinearities in region (1) than in region (2). With respect to these unstable conditions, the FFT spectrum also reflects complex nonlinearities, as observed in www.nature.com/scientificreports/ Fig. 9c. The FFT results from the HBM show the fundamental and relevant subharmonic terms, whereas those of the NS include various harmonic terms. The nonlinear dynamic behavior is more predictable in region (2) than in region (1), but its responses still show that their motions are affected by the period-doubling conditions, even though the HBM is determined to be stable, as shown in Fig. 8. Corresponding to this, the time histories from both the NS and HBM are nearly correlated, as shown in Fig. 9b and d. However, the FFT results of the NS include various harmonic terms around at ̟ 3 (= 46.6 Hz) because its region is still affected by bifurcation, as shown in the marked area (2) of Fig. 8. The dynamic characteristics of the phase diagrams and Poincaré maps are shown in Fig. 10. Figure 10a and c show the dynamic behaviors of region (1), which are significantly affected by complex responses, as shown in Fig. 8. Meanwhile, Fig. 10b and d effectively reflect the dynamic motions in region (2) with one line of the track; however, the Poincare map shows that the periodic motions are not concentrated at one point because their responses are affected by the period-doubling bifurcation, as seen in Fig. 8.  www.nature.com/scientificreports/ Based on the nonlinear frequency response analysis, we propose the desired frequency response function (gray dashed line) tuned by adaptively changing the bistable stiffness, as shown in Fig. 11. Before resonance, the nonlinear stiffness for Case 2 was used to increase the magnitude (i.e., higher sensitivity) and avoid the superharmonic response in Case 1 (or Case 3). After resonance, the nonlinear stiffness was adaptively shifted into Case 1 (or Case 3) to prevent the sub-harmonic response in Case 1. This nonlinear stiffness was recovered in Case 2 to increase the magnitude ratio (i.e., higher sensitivity) beyond ω > 1.5 . Consequently, the frequency response function can be nearly uniformed after resonance (i.e., broadband) 26,27 . Once we fabricate a prototype that mimics the adaptive bistable stiffness of the hair cell bundle, this new adaptive structure enables the development of metamaterial-like broadband vibration applications, such as vibration-based energy harvesters, despite some technical issues that need to be further examined.

Conclusions
In this study, the dynamic behavior of bistable oscillator with adaptive bistable stiffness, which is similar to the hair cell bundle structure in a frog cochlea, was successfully investigated. The main contributions of this study are summarized as follows.
• First, the bistable nonlinearity using the hyperbolic tangent function was effective to avoid instability induced when the discontinuity is connected with the piecewise linear functions when it combines with the harmonic balance method. • Secondly, we investigated all possible nonlinear dynamic characteristics of bistable oscillator by examining the nonlinear frequency responses, phase diagrams, and Poincare maps. • Lastly, we report for the first time that it is possible to achieve a new means of producing broadband (uniform, flat) frequency response functions by mimicking the adaptive bistable stiffness of hair cell bundles.
Overall, it is necessary to implement the proposed bistable oscillator, primarily focusing with adaptive stiffness switching mechanism although our study provides the initial information through the simulation for designing broadband vibration applications such as vibration-based energy harvesters. The adaptive bistable stiffness will be explored by combining different preload adjustment mechanisms with smart-material-based actuators. For the future research direction, we will continue to address some ongoing issues. In particular, further studies such as in-depth bifurcation analysis by mapping the Floquet multipliers are required to complete the nonlinear frequency-response analysis. www.nature.com/scientificreports/

Data availability
The datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.